Take-home_Ex03

Author

Yuheng Liang

Published

October 28, 2024

Modified

October 28, 2024

Take home exe 03

Data

Importing data and package

package

  • spNetwork, which provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.
  • sf package provides functions to manage, processing, and manipulate Simple Features, a formal geospatial data standard that specifies a storage and access model of spatial geometries such as points, lines, and polygons.
  • tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.
  • spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.
  • raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.
  • maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.
pacman::p_load(sf, raster, spatstat, tmap, tidyverse,sparr,spNetwork,dplyr,animation,stringr)

data

drug_case <- read_csv("data/2016-01-01-2024-06-30-Philippines.csv")%>%
    st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
    st_transform(crs = 32651) %>%
  mutate(event_date = dmy(event_date)) %>%
  mutate(event_month = year*100 + month(event_date)) %>%
  mutate(event_quarter = year*10 + quarter(event_date))%>%
  mutate(quarter = as.numeric(str_sub(event_quarter, 5, 5)))
Rows: 6921 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (23): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl  (8): year, time_precision, iso, latitude, longitude, geo_precision, fat...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

import the boundary data

ph_sf = st_read(dsn = "data/phl_adm_psa_namria_20231106_shp",layer = "phl_admbnda_adm2_psa_namria_20231106")
Reading layer `phl_admbnda_adm2_psa_namria_20231106' from data source 
  `/Users/liangyuhang/Downloads/Maaaaaaaaaark/IS415_g/Take-home_Ex/Take-home_Ex03/data/phl_adm_psa_namria_20231106_shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 88 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 114.2779 ymin: 4.587294 xmax: 126.605 ymax: 21.12189
Geodetic CRS:  WGS 84

chack the crs

st_crs(ph_sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

transform the crs

ph_sf = st_transform(ph_sf,crs = 32651)
st_crs(ph_sf)
Coordinate Reference System:
  User input: EPSG:32651 
  wkt:
PROJCRS["WGS 84 / UTM zone 51N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 51N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 120°E and 126°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Japan. North Korea. Philippines. Russian Federation. South Korea. Taiwan."],
        BBOX[0,120,84,126]],
    ID["EPSG",32651]]
ph_sf_owin <-as.owin(ph_sf)

show the point

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(drug_case)+tm_dots()
tmap_mode('plot')
tmap mode set to plotting
hist(drug_case$year)

hist(drug_case$event_month)

drug_case$fatalities <- as.numeric(drug_case$fatalities)

fatalities_by_year <- drug_case %>%
  group_by(year) %>%
  summarise(
    total_fatalities = sum(fatalities, na.rm = TRUE),
    average_fatalities = mean(fatalities, na.rm = TRUE),
    .groups = "drop"
  )

print(fatalities_by_year)
Simple feature collection with 9 features and 3 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -128531.7 ymin: 518340.4 xmax: 864960.5 ymax: 2055587
Projected CRS: WGS 84 / UTM zone 51N
# A tibble: 9 × 4
   year total_fatalities average_fatalities                             geometry
  <dbl>            <dbl>              <dbl>                     <MULTIPOINT [m]>
1  2016             3308              1.41  ((164353.5 1791964), (175482 169682…
2  2017             1959              1.34  ((31818.56 1079518), (168014.2 1814…
3  2018             1211              1.27  ((-128531.7 834249.2), (31818.56 10…
4  2019              805              1.19  ((144968.9 561932.4), (183707.4 167…
5  2020              657              1.19  ((139808.7 560142.4), (155098.6 563…
6  2021              489              1.23  ((100286.8 518340.4), (100576.5 541…
7  2022              232              0.975 ((142101.2 556765.2), (205861.5 568…
8  2023              212              1.08  ((143374.1 558209.3), (196410.6 165…
9  2024               89              1.02  ((144968.9 561932.4), (163956.8 546…
cases_by_month <- drug_case %>%
  group_by(event_month) %>%
  summarise(case_count = n()) %>%
  arrange(event_month)

# 打印结果
print(cases_by_month)
Simple feature collection with 102 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -128531.7 ymin: 518340.4 xmax: 864960.5 ymax: 2055587
Projected CRS: WGS 84 / UTM zone 51N
# A tibble: 102 × 3
   event_month case_count                                               geometry
         <dbl>      <int>                                       <MULTIPOINT [m]>
 1      201601          8 ((240301.8 1676349), (276759.3 1599498), (280878.8 16…
 2      201602         10 ((271691.3 1636477), (280878.8 1622425), (289627.9 16…
 3      201603          5 ((282467.1 1675069), (290083.6 1620386), (296805.5 16…
 4      201604          5 ((214648.2 1775539), (281131.2 1620542), (289422.4 15…
 5      201605         15 ((234185.5 1602344), (262930.6 1586810), (276428.7 16…
 6      201606         63 ((228748.2 1759518), (229359.3 1947400), (237707.6 19…
 7      201607        640 ((175482 1696824), (177128 1788571), (184974.8 166730…
 8      201608        446 ((164353.5 1791964), (198408.4 1761755), (201594 1648…
 9      201609        416 ((198505.8 1644644), (200467.6 1647355), (202671.1 16…
10      201610        291 ((177128 1788571), (207758.6 1774433), (214648.2 1775…
# ℹ 92 more rows
cases_by_quarter <- drug_case %>%
  group_by(event_quarter) %>%
  summarise(case_count = n()) %>%
  arrange(event_quarter)

print(cases_by_quarter)
Simple feature collection with 34 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -128531.7 ymin: 518340.4 xmax: 864960.5 ymax: 2055587
Projected CRS: WGS 84 / UTM zone 51N
# A tibble: 34 × 3
   event_quarter case_count                                             geometry
           <dbl>      <int>                                     <MULTIPOINT [m]>
 1         20161         23 ((240301.8 1676349), (271691.3 1636477), (276759.3 …
 2         20162         83 ((214648.2 1775539), (228748.2 1759518), (229359.3 …
 3         20163       1502 ((164353.5 1791964), (175482 1696824), (177128 1788…
 4         20164        743 ((177128 1788571), (203171.5 1644808), (207313.1 16…
 5         20171        351 ((175482 1696824), (180421.6 1692547), (207313.1 16…
 6         20172        525 ((180421.6 1692547), (207313.1 1641756), (207758.6 …
 7         20173        474 ((31818.56 1079518), (168014.2 1814430), (176145.2 …
 8         20174        117 ((213492.7 1749719), (228983.3 1890994), (240301.8 …
 9         20181        174 ((213706.6 1838983), (218646.9 1807623), (219031.9 …
10         20182        253 ((-128531.7 834249.2), (97838.28 1143181), (193041.…
# ℹ 24 more rows
ggplot(cases_by_month, aes(x = event_month, y = case_count)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  labs(title = "case_month", x = "month", y = "num") +
  theme_minimal()

ggplot(cases_by_quarter, aes(x = event_quarter, y = case_count)) +
  geom_bar(stat = "identity", fill = "lightgreen") +
  labs(title = "case by q", x = "q", y = "num") +
  theme_minimal()

Excluding the third and fourth quarters of 2016 and the first two or three quarters of 2017 and the first two or two quarters of 2021 the distribution of the other years is fairly even not described as seasonal ## Spatial Point Analysis by year Screening of drug use data for each year

Spatio-Temporal Point Pattern Analysis by year

drug_case_2016 <- drug_case %>%
  filter(year== 2016)
drug_case_2017 <- drug_case %>%
  filter(year== 2017)
drug_case_2018 <- drug_case %>%
  filter(year== 2018)
drug_case_2019 <- drug_case %>%
  filter(year== 2019)
drug_case_2020 <- drug_case %>%
  filter(year== 2020)
drug_case_2021 <- drug_case %>%
  filter(year== 2021)
drug_case_2022 <- drug_case %>%
  filter(year== 2022)
drug_case_2023 <- drug_case %>%
  filter(year== 2023)
drug_case_2024 <- drug_case %>%
  filter(year== 2024)

transform them into ppp objects.

drug_case_2016_ppp <- as.ppp(drug_case_2016)
Warning in as.ppp.sf(drug_case_2016): only first attribute column is used for
marks
drug_case_2017_ppp <- as.ppp(drug_case_2017)
Warning in as.ppp.sf(drug_case_2017): only first attribute column is used for
marks
drug_case_2018_ppp <- as.ppp(drug_case_2018)
Warning in as.ppp.sf(drug_case_2018): only first attribute column is used for
marks
drug_case_2019_ppp <- as.ppp(drug_case_2019)
Warning in as.ppp.sf(drug_case_2019): only first attribute column is used for
marks
drug_case_2020_ppp <- as.ppp(drug_case_2020)
Warning in as.ppp.sf(drug_case_2020): only first attribute column is used for
marks
drug_case_2021_ppp <- as.ppp(drug_case_2021)
Warning in as.ppp.sf(drug_case_2021): only first attribute column is used for
marks
drug_case_2022_ppp <- as.ppp(drug_case_2022)
Warning in as.ppp.sf(drug_case_2022): only first attribute column is used for
marks
drug_case_2023_ppp <- as.ppp(drug_case_2023)
Warning in as.ppp.sf(drug_case_2023): only first attribute column is used for
marks
drug_case_2024_ppp <- as.ppp(drug_case_2024)
Warning in as.ppp.sf(drug_case_2024): only first attribute column is used for
marks
drug_case_2016_ppp <- drug_case_2016_ppp[ph_sf_owin]
drug_case_2017_ppp <- drug_case_2017_ppp[ph_sf_owin]
drug_case_2018_ppp <- drug_case_2018_ppp[ph_sf_owin]
drug_case_2019_ppp <- drug_case_2019_ppp[ph_sf_owin]
drug_case_2020_ppp <- drug_case_2020_ppp[ph_sf_owin]
drug_case_2021_ppp <- drug_case_2021_ppp[ph_sf_owin]
drug_case_2022_ppp <- drug_case_2022_ppp[ph_sf_owin]
drug_case_2023_ppp <- drug_case_2023_ppp[ph_sf_owin]
drug_case_2024_ppp <- drug_case_2024_ppp[ph_sf_owin]

re-scale them to km

drug_case_2016_ppp.km <- rescale.ppp(drug_case_2016_ppp, 1000, "km")
drug_case_2017_ppp.km <- rescale.ppp(drug_case_2017_ppp, 1000, "km")
drug_case_2018_ppp.km <- rescale.ppp(drug_case_2018_ppp, 1000, "km")
drug_case_2019_ppp.km <- rescale.ppp(drug_case_2019_ppp, 1000, "km")
drug_case_2020_ppp.km <- rescale.ppp(drug_case_2020_ppp, 1000, "km")
drug_case_2021_ppp.km <- rescale.ppp(drug_case_2021_ppp, 1000, "km")
drug_case_2022_ppp.km <- rescale.ppp(drug_case_2022_ppp, 1000, "km")
drug_case_2023_ppp.km <- rescale.ppp(drug_case_2023_ppp, 1000, "km")
drug_case_2024_ppp.km <- rescale.ppp(drug_case_2024_ppp, 1000, "km")

drug_case_2016_ppp_owin.csr <- envelope(drug_case_2016_ppp.km, Gest, nsim = 99) drug_case_2017_ppp_owin.csr <- envelope(drug_case_2017_ppp.km, Gest, nsim = 99) drug_case_2018_ppp_owin.csr <- envelope(drug_case_2018_ppp.km, Gest, nsim = 99) drug_case_2019_ppp_owin.csr <- envelope(drug_case_2019_ppp.km, Gest, nsim = 99) drug_case_2020_ppp_owin.csr <- envelope(drug_case_2020_ppp.km, Gest, nsim = 99) drug_case_2021_ppp_owin.csr <- envelope(drug_case_2021_ppp.km, Gest, nsim = 99) drug_case_2022_ppp_owin.csr <- envelope(drug_case_2022_ppp.km, Gest, nsim = 99) drug_case_2023_ppp_owin.csr <- envelope(drug_case_2023_ppp.km, Gest, nsim = 99) drug_case_2024_ppp_owin.csr <- envelope(drug_case_2024_ppp.km, Gest, nsim = 99)

plot(drug_case_2016_ppp_owin.csr) plot(drug_case_2017_ppp_owin.csr) plot(drug_case_2018_ppp_owin.csr) plot(drug_case_2019_ppp_owin.csr) plot(drug_case_2020_ppp_owin.csr) plot(drug_case_2021_ppp_owin.csr) plot(drug_case_2022_ppp_owin.csr) plot(drug_case_2023_ppp_owin.csr) plot(drug_case_2024_ppp_owin.csr)

Spatio-Temporal Point Pattern Analysis by quarter

drug_case_2016_q1 <- drug_case %>%
  filter(year== 2016,event_quarter == 20161)
drug_case_2016_q2 <- drug_case %>%
  filter(year== 2016,event_quarter == 20162)
drug_case_2016_q3 <- drug_case %>%
  filter(year== 2016,event_quarter == 20163)
drug_case_2016_q4 <- drug_case %>%
  filter(year== 2016,event_quarter == 20164)
drug_case_2017_q1 <- drug_case %>%
  filter(year== 2017,event_quarter == 20171)
drug_case_2017_q2 <- drug_case %>%
  filter(year== 2017,event_quarter == 20172)
drug_case_2017_q3 <- drug_case %>%
  filter(year== 2017,event_quarter == 20173)
drug_case_2017_q4 <- drug_case %>%
  filter(year== 2017,event_quarter == 20174)
drug_case_2018_q1 <- drug_case %>%
  filter(year== 2018,event_quarter == 20181)
drug_case_2018_q2 <- drug_case %>%
  filter(year== 2018,event_quarter == 20182)
drug_case_2018_q3 <- drug_case %>%
  filter(year== 2018,event_quarter == 20183)
drug_case_2018_q4 <- drug_case %>%
  filter(year== 2018,event_quarter == 20184)
drug_case_2019_q1 <- drug_case %>%
  filter(year== 2019,event_quarter == 20191)
drug_case_2019_q2 <- drug_case %>%
  filter(year== 2019,event_quarter == 20192)
drug_case_2019_q3 <- drug_case %>%
  filter(year== 2019,event_quarter == 20193)
drug_case_2019_q4 <- drug_case %>%
  filter(year== 2019,event_quarter == 20194)
drug_case_2020_q1 <- drug_case %>%
  filter(year== 2020,event_quarter == 20201)
drug_case_2020_q2 <- drug_case %>%
  filter(year== 2020,event_quarter == 20202)
drug_case_2020_q3 <- drug_case %>%
  filter(year== 2020,event_quarter == 20203)
drug_case_2020_q4 <- drug_case %>%
  filter(year== 2020,event_quarter == 20204)
drug_case_2021_q1 <- drug_case %>%
  filter(year== 2021,event_quarter == 20211)
drug_case_2021_q2 <- drug_case %>%
  filter(year== 2021,event_quarter == 20212)
drug_case_2021_q3 <- drug_case %>%
  filter(year== 2021,event_quarter == 20213)
drug_case_2021_q4 <- drug_case %>%
  filter(year== 2021,event_quarter == 20214)
drug_case_2022_q1 <- drug_case%>%
  filter(year== 2022,event_quarter == 20221)
drug_case_2022_q2 <- drug_case%>%
  filter(year== 2022,event_quarter == 20222)
drug_case_2022_q3 <- drug_case%>%
  filter(year== 2022,event_quarter == 20223)
drug_case_2022_q4 <- drug_case%>%
  filter(year== 2022,event_quarter == 20224)
drug_case_2023_q1 <- drug_case%>%
  filter(year== 2023,event_quarter == 20231)
drug_case_2023_q2 <- drug_case%>%
  filter(year== 2023,event_quarter == 20232)
drug_case_2023_q3 <- drug_case%>%
  filter(year== 2023,event_quarter == 20233)
drug_case_2023_q4 <- drug_case%>%
  filter(year== 2023,event_quarter == 20234)
drug_case_2024_q1 <- drug_case%>%
  filter(year== 2024,event_quarter == 20241)
drug_case_2024_q2 <- drug_case%>%
  filter(year== 2024,event_quarter == 20242)
drug_case_2016_ppp1 <- as.ppp(drug_case_2016_q1)
Warning in as.ppp.sf(drug_case_2016_q1): only first attribute column is used
for marks
drug_case_2016_ppp2 <- as.ppp(drug_case_2016_q2)
Warning in as.ppp.sf(drug_case_2016_q2): only first attribute column is used
for marks
drug_case_2016_ppp3 <- as.ppp(drug_case_2016_q3)
Warning in as.ppp.sf(drug_case_2016_q3): only first attribute column is used
for marks
drug_case_2016_ppp4 <- as.ppp(drug_case_2016_q4)
Warning in as.ppp.sf(drug_case_2016_q4): only first attribute column is used
for marks
drug_case_2017_ppp1 <- as.ppp(drug_case_2017_q1)
Warning in as.ppp.sf(drug_case_2017_q1): only first attribute column is used
for marks
drug_case_2017_ppp2 <- as.ppp(drug_case_2017_q2)
Warning in as.ppp.sf(drug_case_2017_q2): only first attribute column is used
for marks
drug_case_2017_ppp3 <- as.ppp(drug_case_2017_q3)
Warning in as.ppp.sf(drug_case_2017_q3): only first attribute column is used
for marks
drug_case_2017_ppp4 <- as.ppp(drug_case_2017_q4)
Warning in as.ppp.sf(drug_case_2017_q4): only first attribute column is used
for marks
drug_case_2018_ppp1 <- as.ppp(drug_case_2018_q1)
Warning in as.ppp.sf(drug_case_2018_q1): only first attribute column is used
for marks
drug_case_2018_ppp2 <- as.ppp(drug_case_2018_q2)
Warning in as.ppp.sf(drug_case_2018_q2): only first attribute column is used
for marks
drug_case_2018_ppp3 <- as.ppp(drug_case_2018_q3)
Warning in as.ppp.sf(drug_case_2018_q3): only first attribute column is used
for marks
drug_case_2018_ppp4 <- as.ppp(drug_case_2018_q4)
Warning in as.ppp.sf(drug_case_2018_q4): only first attribute column is used
for marks
drug_case_2019_ppp1 <- as.ppp(drug_case_2019_q1)
Warning in as.ppp.sf(drug_case_2019_q1): only first attribute column is used
for marks
drug_case_2019_ppp2 <- as.ppp(drug_case_2019_q2)
Warning in as.ppp.sf(drug_case_2019_q2): only first attribute column is used
for marks
drug_case_2019_ppp3 <- as.ppp(drug_case_2019_q3)
Warning in as.ppp.sf(drug_case_2019_q3): only first attribute column is used
for marks
drug_case_2019_ppp4 <- as.ppp(drug_case_2019_q4)
Warning in as.ppp.sf(drug_case_2019_q4): only first attribute column is used
for marks
drug_case_2020_ppp1 <- as.ppp(drug_case_2020_q1)
Warning in as.ppp.sf(drug_case_2020_q1): only first attribute column is used
for marks
drug_case_2020_ppp2 <- as.ppp(drug_case_2020_q2)
Warning in as.ppp.sf(drug_case_2020_q2): only first attribute column is used
for marks
drug_case_2020_ppp3 <- as.ppp(drug_case_2020_q3)
Warning in as.ppp.sf(drug_case_2020_q3): only first attribute column is used
for marks
drug_case_2020_ppp4 <- as.ppp(drug_case_2020_q4)
Warning in as.ppp.sf(drug_case_2020_q4): only first attribute column is used
for marks
drug_case_2021_ppp1 <- as.ppp(drug_case_2021_q1)
Warning in as.ppp.sf(drug_case_2021_q1): only first attribute column is used
for marks
drug_case_2021_ppp2 <- as.ppp(drug_case_2021_q2)
Warning in as.ppp.sf(drug_case_2021_q2): only first attribute column is used
for marks
drug_case_2021_ppp3 <- as.ppp(drug_case_2021_q3)
Warning in as.ppp.sf(drug_case_2021_q3): only first attribute column is used
for marks
drug_case_2021_ppp4 <- as.ppp(drug_case_2021_q4)
Warning in as.ppp.sf(drug_case_2021_q4): only first attribute column is used
for marks
drug_case_2022_ppp1 <- as.ppp(drug_case_2022_q1)
Warning in as.ppp.sf(drug_case_2022_q1): only first attribute column is used
for marks
drug_case_2022_ppp2 <- as.ppp(drug_case_2022_q2)
Warning in as.ppp.sf(drug_case_2022_q2): only first attribute column is used
for marks
drug_case_2022_ppp3 <- as.ppp(drug_case_2022_q3)
Warning in as.ppp.sf(drug_case_2022_q3): only first attribute column is used
for marks
drug_case_2022_ppp4 <- as.ppp(drug_case_2022_q4)
Warning in as.ppp.sf(drug_case_2022_q4): only first attribute column is used
for marks
drug_case_2023_ppp1 <- as.ppp(drug_case_2023_q1)
Warning in as.ppp.sf(drug_case_2023_q1): only first attribute column is used
for marks
drug_case_2023_ppp2 <- as.ppp(drug_case_2023_q2)
Warning in as.ppp.sf(drug_case_2023_q2): only first attribute column is used
for marks
drug_case_2023_ppp3 <- as.ppp(drug_case_2023_q3)
Warning in as.ppp.sf(drug_case_2023_q3): only first attribute column is used
for marks
drug_case_2023_ppp4 <- as.ppp(drug_case_2023_q4)
Warning in as.ppp.sf(drug_case_2023_q4): only first attribute column is used
for marks
drug_case_2024_ppp1 <- as.ppp(drug_case_2024_q1)
Warning in as.ppp.sf(drug_case_2024_q1): only first attribute column is used
for marks
drug_case_2024_ppp2 <- as.ppp(drug_case_2024_q2)
Warning in as.ppp.sf(drug_case_2024_q2): only first attribute column is used
for marks
drug_case_2016_ppp1_owin <- drug_case_2016_ppp1[ph_sf_owin]
drug_case_2016_ppp2_owin <- drug_case_2016_ppp2[ph_sf_owin]
drug_case_2016_ppp3_owin <- drug_case_2016_ppp3[ph_sf_owin]
drug_case_2016_ppp4_owin <- drug_case_2016_ppp4[ph_sf_owin]
drug_case_2017_ppp1_owin <- drug_case_2017_ppp1[ph_sf_owin]
drug_case_2017_ppp2_owin <- drug_case_2017_ppp2[ph_sf_owin]
drug_case_2017_ppp3_owin <- drug_case_2017_ppp3[ph_sf_owin]
drug_case_2017_ppp4_owin <- drug_case_2017_ppp4[ph_sf_owin]
drug_case_2018_ppp1_owin <- drug_case_2018_ppp1[ph_sf_owin]
drug_case_2018_ppp2_owin <- drug_case_2018_ppp2[ph_sf_owin]
drug_case_2018_ppp3_owin <- drug_case_2018_ppp3[ph_sf_owin]
drug_case_2018_ppp4_owin <- drug_case_2018_ppp4[ph_sf_owin]
drug_case_2019_ppp1_owin <- drug_case_2019_ppp1[ph_sf_owin]
drug_case_2019_ppp2_owin <- drug_case_2019_ppp2[ph_sf_owin]
drug_case_2019_ppp3_owin <- drug_case_2019_ppp3[ph_sf_owin]
drug_case_2019_ppp4_owin <- drug_case_2019_ppp4[ph_sf_owin]
drug_case_2020_ppp1_owin <- drug_case_2020_ppp1[ph_sf_owin]
drug_case_2020_ppp2_owin <- drug_case_2020_ppp2[ph_sf_owin]
drug_case_2020_ppp3_owin <- drug_case_2020_ppp3[ph_sf_owin]
drug_case_2020_ppp4_owin <- drug_case_2020_ppp4[ph_sf_owin]
drug_case_2021_ppp1_owin <- drug_case_2021_ppp1[ph_sf_owin]
drug_case_2021_ppp2_owin <- drug_case_2021_ppp2[ph_sf_owin]
drug_case_2021_ppp3_owin <- drug_case_2021_ppp3[ph_sf_owin]
drug_case_2021_ppp4_owin <- drug_case_2021_ppp4[ph_sf_owin]
drug_case_2022_ppp1_owin <- drug_case_2022_ppp1[ph_sf_owin]
drug_case_2022_ppp2_owin <- drug_case_2022_ppp2[ph_sf_owin]
drug_case_2022_ppp3_owin <- drug_case_2022_ppp3[ph_sf_owin]
drug_case_2022_ppp4_owin <- drug_case_2022_ppp4[ph_sf_owin]
drug_case_2023_ppp1_owin <- drug_case_2023_ppp1[ph_sf_owin]
drug_case_2023_ppp2_owin <- drug_case_2023_ppp2[ph_sf_owin]
drug_case_2023_ppp3_owin <- drug_case_2023_ppp3[ph_sf_owin]
drug_case_2023_ppp4_owin <- drug_case_2023_ppp4[ph_sf_owin]
drug_case_2024_ppp1_owin <- drug_case_2024_ppp1[ph_sf_owin]
drug_case_2024_ppp2_owin <- drug_case_2024_ppp2[ph_sf_owin]
drug_case_2016_ppp1.km <- rescale.ppp(drug_case_2016_ppp1_owin, 1000, "km")
drug_case_2016_ppp2.km <- rescale.ppp(drug_case_2016_ppp2_owin, 1000, "km")
drug_case_2016_ppp3.km <- rescale.ppp(drug_case_2016_ppp3_owin, 1000, "km")
drug_case_2016_ppp4.km <- rescale.ppp(drug_case_2016_ppp4_owin, 1000, "km")
drug_case_2017_ppp1.km <- rescale.ppp(drug_case_2017_ppp1_owin, 1000, "km")
drug_case_2017_ppp2.km <- rescale.ppp(drug_case_2017_ppp2_owin, 1000, "km")
drug_case_2017_ppp3.km <- rescale.ppp(drug_case_2017_ppp3_owin, 1000, "km")
drug_case_2017_ppp4.km <- rescale.ppp(drug_case_2017_ppp4_owin, 1000, "km")
drug_case_2018_ppp1.km <- rescale.ppp(drug_case_2018_ppp1_owin, 1000, "km")
drug_case_2018_ppp2.km <- rescale.ppp(drug_case_2018_ppp2_owin, 1000, "km")
drug_case_2018_ppp3.km <- rescale.ppp(drug_case_2018_ppp3_owin, 1000, "km")
drug_case_2018_ppp4.km <- rescale.ppp(drug_case_2018_ppp4_owin, 1000, "km")
drug_case_2019_ppp1.km <- rescale.ppp(drug_case_2019_ppp1_owin, 1000, "km")
drug_case_2019_ppp2.km <- rescale.ppp(drug_case_2019_ppp2_owin, 1000, "km")
drug_case_2019_ppp3.km <- rescale.ppp(drug_case_2019_ppp3_owin, 1000, "km")
drug_case_2019_ppp4.km <- rescale.ppp(drug_case_2019_ppp4_owin, 1000, "km")
drug_case_2020_ppp1.km <- rescale.ppp(drug_case_2020_ppp1_owin, 1000, "km")
drug_case_2020_ppp2.km <- rescale.ppp(drug_case_2020_ppp2_owin, 1000, "km")
drug_case_2020_ppp3.km <- rescale.ppp(drug_case_2020_ppp3_owin, 1000, "km")
drug_case_2020_ppp4.km <- rescale.ppp(drug_case_2020_ppp4_owin, 1000, "km")
drug_case_2021_ppp1.km <- rescale.ppp(drug_case_2021_ppp1_owin, 1000, "km")
drug_case_2021_ppp2.km <- rescale.ppp(drug_case_2021_ppp2_owin, 1000, "km")
drug_case_2021_ppp3.km <- rescale.ppp(drug_case_2021_ppp3_owin, 1000, "km")
drug_case_2021_ppp4.km <- rescale.ppp(drug_case_2021_ppp4_owin, 1000, "km")
drug_case_2022_ppp1.km <- rescale.ppp(drug_case_2022_ppp1_owin, 1000, "km")
drug_case_2022_ppp2.km <- rescale.ppp(drug_case_2022_ppp2_owin, 1000, "km")
drug_case_2022_ppp3.km <- rescale.ppp(drug_case_2022_ppp3_owin, 1000, "km")
drug_case_2022_ppp4.km <- rescale.ppp(drug_case_2022_ppp4_owin, 1000, "km")
drug_case_2023_ppp1.km <- rescale.ppp(drug_case_2023_ppp1_owin, 1000, "km")
drug_case_2023_ppp2.km <- rescale.ppp(drug_case_2023_ppp2_owin, 1000, "km")
drug_case_2023_ppp3.km <- rescale.ppp(drug_case_2023_ppp3_owin, 1000, "km")
drug_case_2023_ppp4.km <- rescale.ppp(drug_case_2023_ppp4_owin, 1000, "km")
drug_case_2024_ppp1.km <- rescale.ppp(drug_case_2024_ppp1_owin, 1000, "km")
drug_case_2024_ppp2.km <- rescale.ppp(drug_case_2024_ppp2_owin, 1000, "km")

drug_case_2016_ppp1.km.csr <- envelope(drug_case_2016_ppp1.km, Gest, nsim = 99) drug_case_2016_ppp2.km.csr <- envelope(drug_case_2016_ppp1.km, Gest, nsim = 99) drug_case_2016_ppp3.km.csr <- envelope(drug_case_2016_ppp1.km, Gest, nsim = 99) drug_case_2016_ppp4.km.csr <- envelope(drug_case_2016_ppp1.km, Gest, nsim = 99) drug_case_2017_ppp1.km.csr <- envelope(drug_case_2017_ppp1.km, Gest, nsim = 99) drug_case_2018_ppp2.km.csr <- envelope(drug_case_2018_ppp1.km, Gest, nsim = 99) drug_case_2019_ppp3.km.csr <- envelope(drug_case_2019_ppp1.km, Gest, nsim = 99) drug_case_2020_ppp4.km.csr <- envelope(drug_case_2020_ppp1.km, Gest, nsim = 99) drug_case_2021_ppp1.km.csr <- envelope(drug_case_2021_ppp1.km, Gest, nsim = 99) drug_case_2022_ppp1.km.csr <- envelope(drug_case_2022_ppp1.km, Gest, nsim = 99) drug_case_2023_ppp1.km.csr <- envelope(drug_case_2023_ppp1.km, Gest, nsim = 99) drug_case_2024_ppp1.km.csr <- envelope(drug_case_2024_ppp1.km, Gest, nsim = 99)

plot(drug_case_2016_ppp1.km.csr)

Spatio-Temporal Point Pattern Analysis by month

drug_case_2016_1 <- drug_case %>%
  filter(year== 2016, event_month == 201601)

transform them into ppp objects.

drug_case_2016_1_ppp <- as.ppp(drug_case_2016_1)
Warning in as.ppp.sf(drug_case_2016_1): only first attribute column is used for
marks

combine

drug_case_2016_1_ppp <- drug_case_2016_1_ppp[ph_sf_owin]
drug_case_2016_1_ppp.km <- rescale.ppp(drug_case_2016_1_ppp, 1000, "km")

drug_case_2016_ppp_owin.csr <- envelope(drug_case_2016_ppp.km, Gest, nsim = 99) plot(drug_case_2016_ppp_owin.csr) Filtering data

Spatio-Temporal Point Pattern Analysis by types

Filtering data

Violence_against_civilians <- drug_case%>%
  filter(event_type == "Violence against civilians")
Strategic_developments <- drug_case%>%
  filter(event_type == "Strategic developments")
Battles <- drug_case%>%
  filter(event_type == "Battles")

Violence against civilians by year

Filtering data

Violence_16 = Violence_against_civilians %>% 
  filter(year == 2016) %>% 
  dplyr::select(quarter)
Violence_17 = Violence_against_civilians %>% 
  filter(year == 2017) %>% 
  dplyr::select(quarter)
Violence_18 = Violence_against_civilians %>% 
  filter(year == 2018) %>% 
  dplyr::select(quarter)
Violence_19 = Violence_against_civilians %>% 
  filter(year == 2019) %>% 
  dplyr::select(quarter)
Violence_20 = Violence_against_civilians %>% 
  filter(year == 2020) %>% 
  dplyr::select(quarter)
Violence_21 = Violence_against_civilians %>% 
  filter(year == 2021) %>% 
  dplyr::select(quarter)
Violence_22 = Violence_against_civilians %>% 
  filter(year == 2022) %>% 
  dplyr::select(quarter)
Violence_23 = Violence_against_civilians %>% 
  filter(year == 2023) %>% 
  dplyr::select(quarter)
Violence_24 = Violence_against_civilians %>% 
  filter(year == 2024) %>% 
  dplyr::select(quarter)
Violence_16_ppp<- as.ppp(Violence_16)
Violence_17_ppp<- as.ppp(Violence_17)
Violence_18_ppp<- as.ppp(Violence_18)
Violence_19_ppp<- as.ppp(Violence_19)
Violence_20_ppp<- as.ppp(Violence_20)
Violence_21_ppp<- as.ppp(Violence_21)
Violence_22_ppp<- as.ppp(Violence_22)
Violence_23_ppp<- as.ppp(Violence_23)
Violence_24_ppp<- as.ppp(Violence_24)
Violence_16_ppp = Violence_16_ppp[ph_sf_owin]
Violence_17_ppp = Violence_17_ppp[ph_sf_owin]
Violence_18_ppp = Violence_18_ppp[ph_sf_owin]
Violence_19_ppp = Violence_19_ppp[ph_sf_owin]
Violence_20_ppp = Violence_20_ppp[ph_sf_owin]
Violence_21_ppp = Violence_21_ppp[ph_sf_owin]
Violence_22_ppp = Violence_22_ppp[ph_sf_owin]
Violence_23_ppp = Violence_23_ppp[ph_sf_owin]
Violence_24_ppp = Violence_24_ppp[ph_sf_owin]

Violence_16_ppp.csr = envelope(Violence_16_ppp, Gest, nsim = 99)

Battles

Filtering data

Battles_16 = Battles %>% 
  filter(year == 2016) %>% 
  dplyr::select(quarter)
Battles_17 = Battles %>% 
  filter(year == 2017) %>% 
  dplyr::select(quarter)
Battles_18 = Battles %>% 
  filter(year == 2018) %>% 
  dplyr::select(quarter)
Battles_19 = Battles %>% 
  filter(year == 2019) %>% 
  dplyr::select(quarter)
Battles_20 = Battles %>% 
  filter(year == 2020) %>% 
  dplyr::select(quarter)
Battles_21 = Battles %>% 
  filter(year == 2021) %>% 
  dplyr::select(quarter)
Battles_22 = Battles %>% 
  filter(year == 2022) %>% 
  dplyr::select(quarter)
Battles_23 = Battles %>% 
  filter(year == 2023) %>% 
  dplyr::select(quarter)
Battles_24 = Battles %>% 
  filter(year == 2024) %>% 
  dplyr::select(quarter)
Battles_16_ppp<- as.ppp(Battles_16)
Battles_17_ppp<- as.ppp(Battles_17)
Battles_18_ppp<- as.ppp(Battles_18)
Battles_19_ppp<- as.ppp(Battles_19)
Battles_20_ppp<- as.ppp(Battles_20)
Battles_21_ppp<- as.ppp(Battles_21)
Battles_22_ppp<- as.ppp(Battles_22)
Battles_23_ppp<- as.ppp(Battles_23)
Battles_24_ppp<- as.ppp(Battles_24)
Battles_16_ppp = Battles_16_ppp[ph_sf_owin]
Battles_17_ppp = Battles_17_ppp[ph_sf_owin]
Battles_18_ppp = Battles_18_ppp[ph_sf_owin]
Battles_19_ppp = Battles_19_ppp[ph_sf_owin]
Battles_20_ppp = Battles_20_ppp[ph_sf_owin]
Battles_21_ppp = Battles_21_ppp[ph_sf_owin]
Battles_22_ppp = Battles_22_ppp[ph_sf_owin]
Battles_23_ppp = Battles_23_ppp[ph_sf_owin]
Battles_24_ppp = Battles_24_ppp[ph_sf_owin]

Battles_16_ppp.csr = envelope(Battles_16_ppp, Gest, nsim = 99)

Strategic developments

Strategic_developments_16 = Strategic_developments %>% 
  filter(year == 2016) %>% 
  dplyr::select(quarter)
Strategic_developments_17 = Strategic_developments %>% 
  filter(year == 2017) %>% 
  dplyr::select(quarter)
Strategic_developments_18 = Strategic_developments %>% 
  filter(year == 2018) %>% 
  dplyr::select(quarter)
Strategic_developments_19 = Strategic_developments %>% 
  filter(year == 2019) %>% 
  dplyr::select(quarter)
Strategic_developments_20 = Strategic_developments %>% 
  filter(year == 2020) %>% 
  dplyr::select(quarter)
Strategic_developments_21 = Strategic_developments %>% 
  filter(year == 2021) %>% 
  dplyr::select(quarter)
Strategic_developments_22 = Strategic_developments %>% 
  filter(year == 2022) %>% 
  dplyr::select(quarter)
Strategic_developments_23 = Strategic_developments %>% 
  filter(year == 2023) %>% 
  dplyr::select(quarter)
Strategic_developments_24 = Strategic_developments %>% 
  filter(year == 2024) %>% 
  dplyr::select(quarter)
Strategic_developments_16_ppp <- as.ppp(Strategic_developments_16)
Strategic_developments_17_ppp <- as.ppp(Strategic_developments_17)
Strategic_developments_18_ppp <- as.ppp(Strategic_developments_18)
Strategic_developments_19_ppp <- as.ppp(Strategic_developments_19)
Strategic_developments_20_ppp <- as.ppp(Strategic_developments_20)
Strategic_developments_21_ppp <- as.ppp(Strategic_developments_21)
Strategic_developments_22_ppp <- as.ppp(Strategic_developments_22)
Strategic_developments_23_ppp <- as.ppp(Strategic_developments_23)
Strategic_developments_24_ppp <- as.ppp(Strategic_developments_24)
Strategic_developments_16_ppp = Strategic_developments_16_ppp[ph_sf_owin]
Strategic_developments_17_ppp = Strategic_developments_17_ppp[ph_sf_owin]
Strategic_developments_18_ppp = Strategic_developments_18_ppp[ph_sf_owin]
Strategic_developments_19_ppp = Strategic_developments_19_ppp[ph_sf_owin]
Strategic_developments_20_ppp = Strategic_developments_20_ppp[ph_sf_owin]
Strategic_developments_21_ppp = Strategic_developments_21_ppp[ph_sf_owin]
Strategic_developments_22_ppp = Strategic_developments_22_ppp[ph_sf_owin]
Strategic_developments_23_ppp = Strategic_developments_23_ppp[ph_sf_owin]
Strategic_developments_24_ppp = Strategic_developments_24_ppp[ph_sf_owin]

Strategic_developments_16_ppp.csr = envelope(Strategic_developments_16_ppp, Gest, nsim = 99)